{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils import column_letter_to_index\n",
    "import pandas as pd\n",
    "import semopy\n",
    "import numpy as np \n",
    "import os\n",
    "\n",
    "# 定义Bootstrap函数\n",
    "def bootstrap_sem(model, data, model_spec,n_bootstrap=2000):\n",
    "\n",
    "    if not os.path.exists('../output/SEM_semopy_bootstrap/'):\n",
    "        os.makedirs('../output/SEM_semopy_bootstrap/')\n",
    "    for ix in range(n_bootstrap):\n",
    "        file_name = '../output/SEM_semopy_bootstrap/'+'{:03d}.txt'.format(ix)\n",
    "        file = open(file_name,'w')\n",
    "        # 生成Bootstrap样本\n",
    "        bootstrap_sample = data.sample(frac=1, replace=True)\n",
    "        \n",
    "        # 重新估计模型\n",
    "        bootstrap_model = semopy.Model(model_spec)\n",
    "        bootstrap_res = bootstrap_model.fit(bootstrap_sample)\n",
    "        \n",
    "        # 提取参数估计结果和拟合指标\n",
    "\n",
    "        file.write('-----------------------bootstrap_model--------------------------\\n')\n",
    "        file.write(str(bootstrap_model.inspect()['Estimate']))\n",
    "        file.write('\\n-----------------------bootstrap_model--------------------------\\n')\n",
    "\n",
    "        file.write('\\n-----------------------bootstrap_res--------------------------\\n')\n",
    "        file.write(str(bootstrap_res))\n",
    "        file.write('\\n-----------------------bootstrap_res--------------------------\\n')\n",
    "        file.close()\n",
    "\n",
    "def SEM_w_sempoy(IV_aph_list, MV_aph_list, DV_aph_list):\n",
    "    \n",
    "    columns_to_extract = {}\n",
    "    for i, IV_aph in enumerate(IV_aph_list):\n",
    "        columns_to_extract[df.columns[IV_aph]] = f\"IV{i+1}\"\n",
    "\n",
    "    for i, MV_aph in enumerate(MV_aph_list):\n",
    "        columns_to_extract[df.columns[MV_aph]] = f\"MV{i+1}\"\n",
    "\n",
    "    for i, DV_aph in enumerate(DV_aph_list):\n",
    "        columns_to_extract[df.columns[DV_aph]] = f\"DV{i+1}\"\n",
    "\n",
    "    processed_df = df[list(columns_to_extract.keys())].rename(columns=columns_to_extract)\n",
    "\n",
    "    model_spec = \"\"\n",
    "    for i in range(len(DV_aph_list)):\n",
    "        model_spec += f\"DV{i+1} ~ \"\n",
    "        for k in range(len(IV_aph_list)):\n",
    "            model_spec += f\"IV{k+1} + \"\n",
    "        for j in range(len(MV_aph_list)):\n",
    "            model_spec += f\"MV{j+1} + \"\n",
    "        model_spec = model_spec[:-3] + \"\\n\"  # Remove extra plus sign\n",
    "\n",
    "    for j in range(len(MV_aph_list)):\n",
    "        model_spec += f\"MV{j+1} ~ \"\n",
    "        for l in range(len(IV_aph_list)):\n",
    "            model_spec += f\"IV{l+1} + \"\n",
    "        model_spec = model_spec[:-3] + \"\\n\"  # Remove extra plus sign\n",
    "   \n",
    "    # 实例化并估计模型\n",
    "    model = semopy.Model(model_spec)\n",
    "    res = model.fit(processed_df)\n",
    "\n",
    "    # 输出模型的估计结果和拟合指标\n",
    "    params = model.inspect()\n",
    "\n",
    "    bootstrap_sem(model, processed_df, model_spec)\n",
    "\n",
    "    fit_stats = None\n",
    "    if hasattr(model, 'last_result'):\n",
    "        fit_stats = model.last_result\n",
    "\n",
    "    return res, params, fit_stats, model_spec\n",
    "\n",
    "\n",
    "if __name__=='__main__':\n",
    "\n",
    "    data_path = '../data/raw_data/334份 按选项序号 汇总变量后.xlsx'\n",
    "    df = pd.read_excel(data_path)\n",
    "\n",
    "   #自变量、中间变量、因变量\n",
    "    IV_list = ['AT','AU','AV'] #自变量X的不同维度，添加几个就加英文逗号隔开\n",
    "    MV_list = ['AY','AZ'] #中间变量\n",
    "    DV_list = ['CH','CI'] #因变量\n",
    " #输出结果在SEM_semopy\n",
    " #看最后一个大表格，两个波兰线的不看\n",
    " #路径系数：看estimate一列（未经归一化的路径系数），计算每个路径的百分比（用该路径数除以全部数的和），百分比即为路径系数\n",
    "    #---------------------------------------------------------------------\n",
    "    IV_index = [column_letter_to_index(IV_) for IV_ in IV_list]\n",
    "    MV_index = [column_letter_to_index(MV_) for MV_ in MV_list]\n",
    "    DV_index = [column_letter_to_index(DV_) for DV_ in DV_list]\n",
    "\n",
    "    IV_index_name = [str(each) for each in IV_index]\n",
    "    MV_index_name = [str(each) for each in MV_index]\n",
    "    DV_index_name = [str(each) for each in DV_index]\n",
    "\n",
    "\n",
    "    file_name = '../output/SEM_semopy/'+'IV_{}_MV_{}_DV_{}.txt'.format('-'.join(IV_index_name),'-'.join(MV_index_name),'-'.join(DV_index_name))\n",
    "    writer = open(file_name,'w',encoding='utf-8')   \n",
    "    writer.write('自变量 {} \\n'.format(str(df.columns[IV_index])))\n",
    "    writer.write('中介变量 {} \\n'.format(str(df.columns[MV_index])))\n",
    "    writer.write('因变量 {} \\n'.format(str(df.columns[DV_index])))  \n",
    "    res,params,fit_stats,model_spec = SEM_w_sempoy(IV_index,MV_index,DV_index)\n",
    "    writer.write('\\n---------------------model_spec---------------------: \\n')\n",
    "    writer.write(str(model_spec))\n",
    "    writer.write('\\n---------------------results---------------------: \\n')\n",
    "    writer.write(str(res))\n",
    "    writer.write('\\n params: \\n')\n",
    "    writer.write(str(params))\n",
    "    writer.write(' \\n')\n",
    "    if fit_stats is not None:\n",
    "        writer.write('\\n---------------------fit_stats---------------------: \\n')\n",
    "        writer.write(str(fit_stats))\n",
    "    writer.write('\\n---------------------chatGPT分析结果---------------------: \\n')\n",
    "    writer.close()\n",
    "\n",
    "    \n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "                \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimate参数均值： [1.16179430e-02 2.21732675e-01 2.32318001e-01 1.85931280e-02\n",
      " 2.72297523e-01 1.90741959e-01 2.04639675e-01 6.63088150e-01\n",
      " 6.59640231e-01 7.21655502e-01 1.56553953e+00 1.18927185e-01\n",
      " 3.05939648e-01 5.54127505e-01 2.66262938e-01 5.26345997e-01\n",
      " 6.84617628e-01 6.75420092e-01 7.07207333e+01 1.03340640e+01]\n",
      "Estimate参数标准差： [0.00660264 0.07655206 0.07070062 0.00722742 0.08136699 0.07269939\n",
      " 0.06514215 0.50725234 0.48781113 0.53047228 0.76717463 0.02638329\n",
      " 0.22913142 0.30362888 0.20239661 0.28057848 0.05905892 0.05998263\n",
      " 5.16080509 0.94160677]\n",
      "Params参数均值： [2.04642000e-01 6.63092500e-01 6.59645500e-01 7.21654500e-01\n",
      " 1.56554250e+00 1.18931000e-01 3.05944000e-01 5.54121000e-01\n",
      " 2.66261000e-01 5.26330000e-01 1.16140000e-02 2.21726500e-01\n",
      " 2.32318000e-01 1.86030000e-02 2.72300500e-01 1.90741500e-01\n",
      " 7.07207235e+01 1.03340785e+01 6.84609000e-01 6.75421000e-01]\n",
      "Params参数标准差： [0.06513313 0.50724579 0.48780819 0.53046918 0.76718856 0.02638663\n",
      " 0.22913328 0.30363453 0.202391   0.28058247 0.00660273 0.07655113\n",
      " 0.0706923  0.00723245 0.08136491 0.07270282 5.16080763 0.94160806\n",
      " 0.05905501 0.05998249]\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import re\n",
    "\n",
    "def analyze_bootstrap_estimates(bootstrap_dir):\n",
    "    # 存储所有bootstrap估计的列表\n",
    "    bootstrap_estimates = []\n",
    "    bootstrap_params = []\n",
    "\n",
    "    # 列出所有的bootstrap文件\n",
    "    bootstrap_files = [file for file in os.listdir(bootstrap_dir) if file.endswith('.txt')]\n",
    "    \n",
    "    # 读取每个bootstrap文件并提取估计值\n",
    "    for file_name in bootstrap_files:\n",
    "        with open(os.path.join(bootstrap_dir, file_name), 'r') as file:\n",
    "            content = file.read()\n",
    "\n",
    "            estimates_section = re.findall(r'-----------------------bootstrap_model--------------------------(.*?)-----------------------bootstrap_model--------------------------', content, re.DOTALL)\n",
    "            if estimates_section:\n",
    "                estimates = re.findall(r'\\b-?\\d+\\.\\d+\\b', estimates_section[0])\n",
    "                estimates = [float(e) for e in estimates]\n",
    "                bootstrap_estimates.append(estimates)\n",
    "\n",
    "            # 使用正则表达式匹配Params的值\n",
    "            params = re.findall(r'(?<=Params: )(.*?)(?=\\n)', content)\n",
    "            params = re.findall(r'\\b-?\\d+\\.\\d+\\b', params[0])\n",
    "            params = [float(e) for e in params]\n",
    "            bootstrap_params.append(params)\n",
    "    \n",
    "    # 将列表转换为NumPy数组以便于计算\n",
    "    bootstrap_estimates = np.array(bootstrap_estimates)\n",
    "    bootstrap_params = np.array(bootstrap_params)\n",
    "    \n",
    "    # 计算每个参数的均值和标准差\n",
    "    estimate_means = np.mean(bootstrap_estimates, axis=0)\n",
    "    estimate_stddevs = np.std(bootstrap_estimates, axis=0)\n",
    "    param_means = np.mean(bootstrap_params, axis=0)\n",
    "    param_stddevs = np.std(bootstrap_params, axis=0)\n",
    "    \n",
    "    return estimate_means, estimate_stddevs, param_means, param_stddevs\n",
    "\n",
    "# 使用实例\n",
    "bootstrap_dir = '../output/SEM_semopy_bootstrap/'\n",
    "estimate_means, estimate_stddevs, param_means, param_stddevs = analyze_bootstrap_estimates(bootstrap_dir)\n",
    "print(\"Estimate参数均值：\", estimate_means)\n",
    "print(\"Estimate参数标准差：\", estimate_stddevs)\n",
    "print(\"Params参数均值：\", param_means)\n",
    "print(\"Params参数标准差：\", param_stddevs)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "whj",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
